Statistik och dataanalys I, 15 hp

Datorlaboration 7 - Inferens för regression

Mattias Villani

Publicerad

October 22, 2023

Installera paket

Den här labben förutsätter att följande paket finns installerade:

  • mosaic
  • sda123

sda123 är kursens egna R-paket och måste installeras från GitHub. Det görs genom att först installera paketet remotes och därefter installera kurspaketet:

library(mosaic)
# install.packages("remotes") # avkommentera första gången
library(remotes)
#install_github("StatisticsSU/sda1paket")  # avkommentera första gången
library(sda123)

Introduktion

I den här datorlabben kommer ni analysera olika datamaterial med regression. Det första datamaterialet om reklam är verkligt. I Uppgift 3 ska ni analysera datamaterial som är simulerade av mig, från olika populationsmodeller som jag har valt. Det blir ett intressant laboratorium, där ni kan utforska olika aspekter av regression i en kontrollerad miljö. Och även få veta den underliggande populationsmodellen efter att ni gjort labben! The truth will reveal itself!

Skapa mapp för labben

Skapa en mapp Lab7 i din kursmapp SDA1. Ladda ner Quarto-filen för denna lab genom att högerklicka här och välj ‘Spara länk’ eller något likande från menyn som dyker upp. Spara filen i din Lab7 mapp. Öppna Quarto-filen i RStudio och fortsätt med denna laboration direkt i Quarto-dokumentet, där du också fyller i svaren på dina laborationsövningar.

1. Enkel regression för reklamdata

Datamaterialet reklam.csv innehåller data på n=200 produkters försäljning (sales) i tusentals US dollar och hur mycket (i s k reklamenheter) produkten har marknadsförts i olika reklamkanaler: tv, radio och newspaper. Vi läser in data med kommandot:

reklam = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/reklam.csv')
head(reklam)
     tv radio newspaper sales
1 230.1  37.8      69.2  22.1
2  44.5  39.3      45.1  10.4
3  17.2  45.9      69.3  12.0
4 151.5  41.3      58.5  16.5
5 180.8  10.8      58.4  17.9
6   8.7  48.9      75.0   7.2

Vi börjar med en enkel regressionsmodell med enbart tv som förklarande variabel. Vi utforskar data genom histogram för responsvariabeln sales och scatterplot mot tv.

par(mfrow = c(1,2))
hist(reklam$sales, col = "orange", main = "")
plot(sales ~ tv, data = reklam, pch = 19, cex = 0.7, col = "steelblue")

💪 Uppgift 1.1

Skatta den enkla linjära regressionsmodellen i R:

\[ \texttt{sales} = \beta_0 + \beta_1 \cdot \texttt{tv} + \varepsilon, \hspace{1cm} \varepsilon\overset{iid}{\sim}N(0,\sigma_{\varepsilon}) \]

och tolka skattningen av intercept och lutningskoefficient. Vilken skattning av \(\sigma_{\varepsilon}\) ger R, och hur tolkar du denna?

Den skattningen av intercept är ca 6.97 och lutningskoefficient är ca 0.06 och skattning av \(\sigma_{\varepsilon}\) är ca 2.30.

model1 = lm(sales ~ tv, data = reklam)
summary(model1)

Call:
lm(formula = sales ~ tv, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4438 -1.4857  0.0218  1.5042  5.6932 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.974821   0.322553   21.62   <2e-16 ***
tv          0.055465   0.001896   29.26   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared:  0.8122,    Adjusted R-squared:  0.8112 
F-statistic: 856.2 on 1 and 198 DF,  p-value: < 2.2e-16

💪 Uppgift 1.2

Testa om variabeln tv är en signifikant förklarande variabel på 5% signifikansnivå. Ställ upp noll- och alternativhypotes, teststatistiska med fördelning under nollhypotesen. Utför testet med hjälp av formelsamlingen och formulera din slutsats. Du får använda information från R, men du ska ställa upp formeln för teststatistikan och beräkna med den formeln.

model1 = lm(sales ~ tv, data = reklam)
confint(model1)
                 2.5 %     97.5 %
(Intercept) 6.33874038 7.61090260
tv          0.05172671 0.05920283
summary(model1, conf_intervals = TRUE, anova = FALSE)

Call:
lm(formula = sales ~ tv, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4438 -1.4857  0.0218  1.5042  5.6932 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.974821   0.322553   21.62   <2e-16 ***
tv          0.055465   0.001896   29.26   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared:  0.8122,    Adjusted R-squared:  0.8112 
F-statistic: 856.2 on 1 and 198 DF,  p-value: < 2.2e-16
t_teststatistikan <- (0.06-0)/0.002
t_teststatistikan
[1] 30
n = nrow(reklam)

qt(0.975,n-2)
[1] 1.972017
t_teststatistikan > qt(0.975,n-2) # if TRUE -> reject H0: Beta1 = 0
[1] TRUE

💪 Uppgift 1.3

Visa hur R har beräknat p-värdet i regression-utskriften (från summary) genom att göra beräkningen själv med relevant s k p,d,q, eller r-funktion.

model1 = lm(sales ~ tv, data = reklam)
summary(model1)

Call:
lm(formula = sales ~ tv, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4438 -1.4857  0.0218  1.5042  5.6932 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.974821   0.322553   21.62   <2e-16 ***
tv          0.055465   0.001896   29.26   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared:  0.8122,    Adjusted R-squared:  0.8112 
F-statistic: 856.2 on 1 and 198 DF,  p-value: < 2.2e-16
p_value <- 2*pt(t_teststatistikan, n-2, lower.tail = F) #P(Tobs>=1.97)
p_value
[1] 1.399538e-75
pt(1.97, df=198)#sannolikhet
[1] 0.9748837

💪 Uppgift 1.4

Gör en residualanalys med funktionen reg_residuals i sda123-paketet, och kommentera om var och ett av de fyra grundläggande antaganden för populationsmodellen verkar uppfyllda.

Figur ett är qqplot, data ligger på theiretical linje och tolkar som normalfördelning.

library(sda123)

model1 = lm(sales ~ tv, data = reklam)
summary(model1)

Call:
lm(formula = sales ~ tv, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4438 -1.4857  0.0218  1.5042  5.6932 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.974821   0.322553   21.62   <2e-16 ***
tv          0.055465   0.001896   29.26   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared:  0.8122,    Adjusted R-squared:  0.8112 
F-statistic: 856.2 on 1 and 198 DF,  p-value: < 2.2e-16
reg_residuals(model1)

💪 Uppgift 1.5

Gör en prediktion med 95%-igt prediktionsintervall för sales vid tv=100. [tips: argumentet newdata i predict-funktionen måste vara en dataframe, inte en vektor. Se min kod för lifespan data.]

model1 = lm(sales ~ tv, data = reklam)
reg_summary(model1)

Analysis of variance - ANOVA
------------------------------------------------
       df     SS        MS      F     Pr(>F)
Regr    1 4512.4 4512.4352 856.18 7.9279e-74
Error 198 1043.5    5.2704                  
Total 199 5556.0                            

Measures of model fit
------------------------------------------------
Root MSE       R2   R2-adj 
 2.29575  0.81218  0.81123 

Parameter estimates
------------------------------------------------
            Estimate Std. Error t value   Pr(>|t|)
(Intercept) 6.974821  0.3225535  21.624 5.0277e-54
tv          0.055465  0.0018956  29.260 7.9279e-74
# Prediktionsintervall - sda123
reg_predict(sales ~ tv, data = reklam)
`geom_smooth()` using formula = 'y ~ x'

qt(p=0.975, df=198) #sannolikhet
[1] 1.972017
# Anpassat värde/prediktion för Sverige:
predict(model1, newdata = data.frame(tv = 100),interval="prediction")
      fit      lwr      upr
1 12.5213 7.979338 17.06326

2. Multipel regression för reklamdata

Vi går nu vidare med en multipel regressionsmodell med alla tre förklarande variabler: tv, radio och newspaper. Först en titt på data:

par(mfrow = c(2,2))
hist(reklam$sales, col = "orange", main = "")
plot(sales ~ tv, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
plot(sales ~ radio, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
plot(sales ~ newspaper, data = reklam, pch = 19, cex = 0.7, col = "steelblue")

💪 Uppgift 2.1

Anpassa följande multipla regressionsmodell:

\[ \texttt{sales} = \beta_0 + \beta_1 \cdot \texttt{tv} + \beta_2 \cdot \texttt{radio} + \beta_3 \cdot \texttt{newspaper} + \varepsilon, \hspace{1cm} \varepsilon\overset{iid}{\sim}N(0,\sigma_{\varepsilon}) \]

och tolka skattningen av regressionskoefficienten för tv.

model2 <- lm(sales ~ tv + radio + newspaper, data = reklam)
summary(model2)

Call:
lm(formula = sales ~ tv + radio + newspaper, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.3034 -0.8244 -0.0008  0.8976  3.7473 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.6251241  0.3075012  15.041   <2e-16 ***
tv          0.0544458  0.0013752  39.592   <2e-16 ***
radio       0.1070012  0.0084896  12.604   <2e-16 ***
newspaper   0.0003357  0.0057881   0.058    0.954    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.662 on 196 degrees of freedom
Multiple R-squared:  0.9026,    Adjusted R-squared:  0.9011 
F-statistic: 605.4 on 3 and 196 DF,  p-value: < 2.2e-16

💪 Uppgift 2.2

Testa om variabeln tv är en signifikant förklarande variabel på 5% signifikansnivå. Ställ upp noll- och alternativhypotes, teststatistiska med fördelning under nollhypotesen. Utför testet och formulera din slutsats. Du får använda information från R, men du ska ställa upp formeln för teststatistikan och beräkna den.

Kommentera kortfattat utifrån utskriften om radio och newspaper är signifikanta på 5% signifikansnivå.

Radio är signifikanta på 5% signifikansnivå utan newspaper inte är signifikanta på 5% signifikansnivå.

model2 <- lm(sales ~ tv + radio + newspaper, data = reklam)
confint(model2)
                  2.5 %     97.5 %
(Intercept)  4.01868836 5.23155980
tv           0.05173372 0.05715784
radio        0.09025861 0.12374384
newspaper   -0.01107921 0.01175052
summary(model2)

Call:
lm(formula = sales ~ tv + radio + newspaper, data = reklam)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.3034 -0.8244 -0.0008  0.8976  3.7473 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.6251241  0.3075012  15.041   <2e-16 ***
tv          0.0544458  0.0013752  39.592   <2e-16 ***
radio       0.1070012  0.0084896  12.604   <2e-16 ***
newspaper   0.0003357  0.0057881   0.058    0.954    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.662 on 196 degrees of freedom
Multiple R-squared:  0.9026,    Adjusted R-squared:  0.9011 
F-statistic: 605.4 on 3 and 196 DF,  p-value: < 2.2e-16
# H0: beta_tv = 0, Ha: beta_tv != 0
t_teststatistikan <- (0.0544458-0)/0.0013752
t_teststatistikan
[1] 39.59119
n = nrow(reklam)

qt(0.975,n-2)
[1] 1.972017
t_teststatistikan > qt(0.975,n-2) # if TRUE -> reject H0: Beta1 = 0
[1] TRUE

💪 Uppgift 2.3

Vilken modell, den enkla regressionen eller den multipla, föredrar du? Motivera med relevant R² värde och genom att göra 5-fold korsvalidering med funktionen reg_crossval i sda123-paketet.

RMSE_model2 = reg_crossval(sales ~ tv + radio + newspaper, data = reklam, nfolds = 5)
RMSE_model2 
[1] 1.677374

3. Residualanalys på simulerade data

Jag har simulerat fem olika datamaterial simdata1-5 från en linjär regressionsmodell:

\[ \texttt{y} = \beta_0 + \beta_1 \cdot \texttt{X1} + \beta_2 \cdot \texttt{X2} + \varepsilon \]

där i vissa datamaterial har feltermer \(\varepsilon\) som inte uppfyller alla av de fyra antagandena, och eventuellt kan det också finnas ett icke-linjärt samband mellan någon förklarande variabel och y. Er uppgift är att utföra residualanalys med funktionen reg_residuals i sda123-paketet för att försöka upptäcka vilka antaganden som inte är uppfyllda i respektive datamaterial (Uppgift 3.1-3.5 nedan). Ni bör också titta på plot(simdata1) för varje datamaterial för att upptäcka icke-linjära samband.

Här är de fem datamaterialen:

simdata1 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata1.csv')
simdata2 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata2.csv')
simdata3 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata3.csv')
simdata4 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata4.csv')
simdata5 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata5.csv')

💪 Uppgift 3.1 - simdata1

head(simdata1)
  X          y          X1          X2
1 1  1.2856459  0.30665435  0.66858877
2 2  0.4728535 -0.59227811  0.29043835
3 3  0.4576525 -1.48905857  0.21114563
4 4 -0.1214482 -0.64102846 -0.55785480
5 5 -3.2053896 -1.20232219 -0.09721979
6 6 -0.3920066  0.04123726  1.71397704
m1 = lm(y ~ X1 + X2, data = simdata1)
reg_residuals(m1)

summary(m1)

Call:
lm(formula = y ~ X1 + X2, data = simdata1)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.3951  -0.5728   0.1052   0.8056   5.6914 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.9851     0.1101   8.945 2.74e-16 ***
X1            1.0369     0.1035  10.014  < 2e-16 ***
X2           -0.1385     0.1067  -1.298    0.196    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.541 on 197 degrees of freedom
Multiple R-squared:  0.3513,    Adjusted R-squared:  0.3447 
F-statistic: 53.35 on 2 and 197 DF,  p-value: < 2.2e-16
plot(simdata1)

💪 Uppgift 3.2 - simdata2

m2 = lm(y ~ X1 + X2, data = simdata2)
reg_residuals(m2)

summary(m2)

Call:
lm(formula = y ~ X1 + X2, data = simdata2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5471 -0.5977  0.0276  0.6443  4.0607 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03782    0.07845  13.230   <2e-16 ***
X1           0.99503    0.08783  11.329   <2e-16 ***
X2           0.04638    0.07719   0.601    0.549    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.108 on 197 degrees of freedom
Multiple R-squared:  0.3952,    Adjusted R-squared:  0.3891 
F-statistic: 64.38 on 2 and 197 DF,  p-value: < 2.2e-16
plot(simdata2)

💪 Uppgift 3.3 - simdata3

m3 = lm(y ~ X1 + X2, data = simdata3)
reg_residuals(m3)

summary(m3)

Call:
lm(formula = y ~ X1 + X2, data = simdata3)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.67225 -0.53425  0.09314  0.55420  2.02347 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03070    0.06272  16.434   <2e-16 ***
X1           0.92495    0.06195  14.931   <2e-16 ***
X2           0.09723    0.06693   1.453    0.148    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8854 on 197 degrees of freedom
Multiple R-squared:  0.5352,    Adjusted R-squared:  0.5305 
F-statistic: 113.4 on 2 and 197 DF,  p-value: < 2.2e-16
plot(simdata3)

💪 Uppgift 3.4 - simdata4

m4 = lm(y ~ X1 + X2, data = simdata4)
reg_residuals(m4)

summary(m4)

Call:
lm(formula = y ~ X1 + X2, data = simdata4)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3459 -1.4186 -0.0414  1.3921  5.0268 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.68543    0.14816   4.626 6.73e-06 ***
X1           0.83993    0.13810   6.082 6.07e-09 ***
X2           0.02523    0.13708   0.184    0.854    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.077 on 197 degrees of freedom
Multiple R-squared:  0.1582,    Adjusted R-squared:  0.1496 
F-statistic: 18.51 on 2 and 197 DF,  p-value: 4.302e-08
plot(simdata4)

💪 Uppgift 3.5 - simdata5

m5 = lm(y ~ X1 + X2, data = simdata5)
reg_residuals(m5)

summary(m5)

Call:
lm(formula = y ~ X1 + X2, data = simdata5)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.8418 -0.2863  0.0024  0.2579  1.0555 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.35010    0.02720  12.873   <2e-16 ***
X1           0.02583    0.04673   0.553    0.581    
X2          -0.03729    0.02755  -1.354    0.177    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3824 on 197 degrees of freedom
Multiple R-squared:  0.01134,   Adjusted R-squared:  0.001303 
F-statistic:  1.13 on 2 and 197 DF,  p-value: 0.3252
plot(simdata5)